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Abstract. We develop a new theory for treating boundary problems 
for linear ordinary differential equations whose fundamental system may 
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1 Introduction 

The treatment of boundary problems in symbolic computation was initiated 
in the PhD thesis [S] under the guidance of Bruno Buchberger in cooperation 
with Heinz Engl; see also [HH] and for the further development mm- Its 
implementation was originally carried out within the TH30REMV project [T]. 

Up to now, we have always assumed differential equations without singularity 
or, equivalently, monic differential operators (leading coefficient function being 
unity). In this paper, we develop for the first time an algebraic theory for treating 
boundary problems with a (mild) singularity at one endpoint. Eor details, we 
refer to Section Our approach is very different from the traditional analysis 
setting in terms of the Weyl-Titchmarsh theory (limit points and limit circles). 
It would be very interesting to explore the connections between our approach 
and the classical treatment; however, this must be left for future work. 

Regarding the general setup of the algebraic language for boundary problems, 
we refer to the references mentioned above, in particular m- At this point, 
let us just recall some notation. We start from a fixed integro-differential al¬ 
gebra (J^, 9, J). The formulation of (local) boundary conditions is in terms of 
evaluations, which are by definition multiplicative linear functionals in T* . We 
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write J-\ < J-2 if is a subspace of J-2; the same notation is used for subspaces 
of the dual space JF* . The orthogonal of a subspace B < J-* is defined as 

B^ = {uGT\ P{u) = 0 for all P e B} 

and similarly for the orthogonal of a subspace A. < A. We write 
for the linear span of (possibly infinitely many) elements /i,/ 2 r’’ G A] the 
same notation is used for linear spans within A* . The zero vector space of any 
dimension is denoted by O. 

We write 

n- = n{n — 1) • • • (n — fc + 1) 

for the falling factorial, where n S C could be arbitrary (but will be an integer for 
our purposes) while k is taken to be a natural number. Note that n- = 0 for k > 
n. Our main example of an integro-differential algebra will be the space C‘^{I) 
of analytic functions on a closed interval / = [a, 6]. Recall that by definition this 
is the space of functions that are analytic in some open set containing [a, b\. 

The development of the new algebraic theory of two-point boundary prob¬ 
lems with a mild singularity (whose treatment is really just broached in this 
paper) was triggered by a collaboration between a symbolic computation team 
(consisting of the first, second and fourth author) and a researcher in engineering 
mechancs (the third author). This underlines the importance and fruitfulness of 
collaborations between theoretical developments and practical applications. We 
present the example that had originally led to this research in Section]^ 

From a methodological point of view, this research on the symbolic algorithm 
for linear boundary value problems has particular relevance and is drawing from 
the Theorema Project, see [T]. This project aims at supporting a new paradigm 
for doing (algorithmic) mathematical research: In the phase of doing research 
on new theorems and algorithms, TH30REMV provides a formal language (a 
version of predicate logic) and an automated reasoning system by which the 
exploration of the theory is supported. In the phase in which algorithms based 
on the new theory should be implemented and used in computing examples, 
TH30REMV allows to program and execute the algorithms in the same lan¬ 
guage. In the particular case of our approach to solving linear boundary value 
problems, the fundamental theorem on which the approach is based was proved 
automatically by checking that the rewrite rules for integro-differential opera¬ 
tors forms a Grobner basis. In a second step, the algorithm for solving linear 
boundary value problems is expressed again in the TH30REMV language and 
can then be called by the users by inputting the linear boundary value problems 
in a user-friendly notation. 

In its current version, the engine for solving boundary problems is bundled 
in the GreenGroebner package of TH30REMV. As an example, consider the 
boundary problem 

u(0) = u(l) = 0, 
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which we shall consider in greater detail in Section This can be given to 
TH30REMV in the form 


BPSolve 


u" -h i u' - i u = f 
u[0] = u[l] = 0 


,u, X, 0, 


1 


leading either to the solution for u{x) as 

or to the Green’s function g{x,^) as 

rii(-i+x2)f ^ e<x 
\ix(-i + e') ^ x<e 

at the user’s request. 

In addition to the GreenGroebner package in TH30REMV, a Maple package 
named IntDiffOp is also available [5]. This package was developed in the frame 
of Anja Korporal’s PhD thesis, supervised by Georg Regensburger and the first 
author. The Maple package supports also generalized boundary problems (see 
Section]^ for their relevance to this paper). One advantage of the TH30REMV 
system is that both the research phase and the application phase of our method 
can be formulated and supported within the same logic and software system— 
which we consider to be quite a novel and promising paradigm for the future. 


2 A Simple Example 


For illustrating the new ideas, it is illuminating to look at a simple example that 
exhibits the kind of phenomena that we have to cope with in the Kirchhoff plate 
boundary problem. 


Example 1. Let us start with the intuitive but mathematically unprecise state¬ 
ment of the following example: Given a “suitable” forcing function / on the unit 
interval I = [0,1], we want to find a “reasonable” solution function u such that 


u” + ^u' - ^u = f, 
u(0) = u(l) = 0. 


( 1 ) 


But note that the differential operator T = \s singular at the 

left boundary point x = 0 of the interval I under consideration. Hence the first 
boundary condition u(0) = 0 should be looked at with some suspicion. And what 
function space are we supposed to consider in the first place? If the ^ and ^ are 
to be taken literally, the space C'°°[0,1] will clearly not do. On the other hand, 
we need functions that are smooth (or at least continuous) at x = 0 to make 
sense of u(0). 
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In the rest of this section, we shall give one possible solution to the dilemma 
outlined above. Of course we could resort to using different function space for u 
and /, and this is in fact the approach one usually takes in Analysis. For our 
present purposes, however, we prefer to keep the simple paradigm of integro- 
differential algebras as outlined in Section but we shall modify it to accom¬ 
modate singularities such as in ^ and Note that these are just poles, so we 
can take J" to be the subring of the field M{I) consisting of complex-valued 
meromorphic functions that are regular at all a: S / except possibly a; = 0. In 
other words, these are functions that have a Laurent expansion at a: = 0 with 
finite principal part, converging in a complex annulus 0 < j^l < p with p > 1. 
In fact, we will only use the real part [—1,1] \ {0} of this annulus. Note that we 
have of course y, G J-. 

The ring iF is an integro-differential algebra over C if we use the standard 
derivation d = -^ and the Rota-Baxter operator 

initialized at the regular point a: = 1. 

We have now ensured that the differential operator of Q has a clear algebraic 
interpretation T G F^[d\. However, the boundary condition u(0) = 0 is still 
dubious. For making it precise, note that the integro-differential algebra {F, d, J) 
has only the second of the two boundary evaluations L,R: —>■ C with L{f ) = 
/(O) and R{f) = /(I) in the usual sense of a total function. So while we can 
interpret the second boundary condition algebraically by Ru = 0, the same 
does not work on the left endpoint. Instead of an evaluation at a: = 0 we shall 
introduce the map 

OO —1 

pp: ^ Onx'^ ^ 

n—N n—N 

that extracts the principal part of a function written in terms of its Laurent 
expansion at x = 0. Here and henceforth we assume such expansions of nonzero 
functions are written with ojv 0. If > 0 the function is regular at x = 0, 
and the above sum is to be understood as zero. Clearly, pp: J" —)■ is a linear 
projector, with the complementary projector reg := Ijr — pp extracting the 
regular part at x = 0. Incidentially, pp and reg are also Rota-Baxter operators of 
weight —1, which play a crucial role in the renormalization theory of perturbative 
quantum field theory [3l Ex. 1.1.10]. 

Finally, we define the functional C: A —> C that extracts the constant 
term oq of a meromorphic function expanded at x = 0. Combining C with the 
monomial multiplication operators, we obtain the coefficient functionals (x") := 
Cx”" (n G Z) that map to a„. In particular, the residue functional is 

given by (x“^) = Cx. 

Note that for functions regular at x = 0, the functional C coincides with 
the evaluation at the left endpoint, L: F —>■ C, f i-G /(O). However, for general 
meromorphic functions, L is undefined and C is not multiplicative since for 
example C(x-y) = 1^0 = C'(x) C'(y). Hence we refer to C only as a functional 
but not as an evaluation. 
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We can now make the boundary condition u(0) =0 precise for our setting 
over J-. What we really mean is that A(w) = 0, where A := pp-t-C is the projector 
that extracts the principal part together with the constant term. Extending the 
algebraic notion of boundary problem to allow for boundary conditions that are 
not functionals, we may thus view 0 as + Id - [A, R]p. In this way 

we have given a precise meaning to the formulation of the boundary problem. 

But how are we to go about its solution in a systematic manner? Let us first 
look at the adhoc standard method way of doing this—determining the general 
homogeneous solution, then add the inhomogeneous solution via variation of con¬ 
stants, and finally adapt the integration constants to accommodate the boundary 
conditions. In our case, one sees immediately that Ker(T) = [a;, A] so that the 
general solution of the homogeneous differential equation is u{x) = cix + 
where ci,C 2 S C are integration constants. Variation of constants [21 p. 74] then 
yields 

u{x) = cix + ^ + (2) 

for the inhomogeneous solution. Note that f G R may also have singularities 
as a: = 0 or any other point x G I apart from x = 1. 

Now we need to impose the boundary conditions. From w(l) = 0 we obtain 
immediately Ci -I-C 2 = 0. For the boundary condition at x = 0 we have to proceed 
a bit more cautiously, obtaining 

„(o) = lim(c.x+1 +1^ nodi-^l enodi) 


where we assume that / is regular at 0. It is clear that the remaining limit can 
only exist if the integral tends to a finite limit as x —)■ 0, and this is the case 
by our assumption on /. But then we may apply the boundary condition to 
obtain C 2 = | df. This gives the overall solution 


u{x) 




1 

2x 


e]fpd^, 


( 3 ) 


which one may write in the standard form u{x) = Jg 5 (x, ^)/(^) where the 
Green’s function is defined as 





if ^ < X 
if ^ > X 


( 4 ) 


in the usual manner. 

How are we to make sense of this in an algebraic way, i.e. without (explicit) 
use of limits and hence topology? The key to this lies in the projector pp and 
the functional L, which serve to distill into our algebraic setting what we need 
from the topology (namely /(x) = (pp/)(x) -I- 0(1) as x —)■ 0). However, there 
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is another complication when compared to boundary problems without singu¬ 
larities as presented in Section]^ We cannot expect a solution u G to Q for 
every given forcing function / G In other words, this boundary problem is 
not regular in the sense of [la¬ 
in the past, we have also used the term singular boundary problem for such 
situations (here this seems to be suitable in a double sense but we shall be careful 
to separate the second sense by sticking to the designation “boundary problems 
with singularities”). The theory of singular boundary problems was developed 
in an abstract setting in |4]; applications to boundary problems (without singu¬ 
larities) have been presented in [^ . At this point we shall only recall a few basic 
facts. 

A boundary problem (T, B) is called semi-regular if Ker(T) H = O. 
It is easy to see that the boundary problem (U^ -|- — ^,[A,i?]-*-) is in 

fact semi-regular. Since any u G Ker(T) can be written as u{x) = cix , 
the condition Ru = 0 implies C 2 = —ci and hence u{x) = ci{x — ^). But 
then (Au)(x) = — ^ = 0 forces ci = 0 and hence u = 0. 

If (ZI^ -I- — ^,[A, i?]-*-) were a regular boundary problem, we would 

have Ker(T) -j- [A, i?]-*- = R. However, it is easy to see that there are ele¬ 
ments u G T that do not belong to Ker(T) -|- [A,!?]-*", for example u{x) = 
Hence we conclude that the boundary problem 0 is in fact overdetermined. For 
such boundary problems (T, B) one can always select a regular subproblem (T, B), 
in the sense that In our case, a natural choice is = [(x“^), i?]. This is 

regular since the evaluation matrix 

((x-i),i?)(i,x) = Q 

is regular, and the associated kernel projector is P = ^ {x~^) x {R — {x~^)) 
by [3 Lem. A.l]. 

For making the boundary problem Q well-defined on P we need one more 
ingredient: We have to fix a complement £ of T{B'^), the so-called exceptional 
space. Intuitively speaking, this comprises the “exceptional functions” of R that 
we decide to discard in order to render 0 solvable. Let us first work out 
what T{B^) < IF looks like. Since u{x) = S B^ forces the princi¬ 

pal part and constant term to vanish, we may start from u{x) = 
with the additional proviso that X]n>o = 1- Applying T — D'^ -\- to 

this u{x) yields — l)x"'“^, which represents an arbitrary element 

of C^{I) = [pp]''" for a suitable choice of coefficients (a„)„>i since the addi¬ 
tional condition J2n>o ^ can always be met by choosing oi = 1 — X]ra>i 
But then it is very natural to choose £ = [reg]-*- as the required complement. 
Clearly, the elements of this space £ are the Laurent polynomials of C[i] without 
constant term. 

By Prop. 2 of [S], the Green’s operator of a generalized boundary prob¬ 
lem {T,B,£) is given by G = GQ, where G is the Green’s operator of some 
regular subproblem (T, B) and Q is the projector onto T{B^) along £. In our case, 
the latter projector is clearly Q = reg, while the Green’s operator G = (I —P)T^ 
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by [111 Thm. 26], with the kernel projector P as above and the fundamental right 
inverse given according to [TTJ Prop. 23] by 

= I Aix - ^ Aix'^, 

which is essentially just a reformulation of the inhomogeneous part in Q. Fol¬ 
lowing the style of [9] we write here as Ai G J] to emphasize its role 
in the integro-differential operator ring. Similarly, we write F := —LAi = 
for the definite integral over the full interval /, which we may regard as a linear 
functional C^{I) —>■ C. 

Putting things together, it remains to compute 

G = {I - P)T^Q = (1 - ^ -b x - xR){\ Aix - ^ A^x"^) reg 
= ( - f i i Fx^ -b f J^x2)reg, 

which may be done by the usual rewrite rules El Tbl. 1] for the operator 
ring J'[9, J], together with the obvious extra rules that on Im(reg) = C^{I) the 
residual (x“^) vanishes and C = (x°) coincides with the evaluation L at zero. 
Using the standard procedure for extracting the Green’s function, one obtains 
exactly Q when restricting the forcing functions to f G C'“(/). We have thus 
succeeded in applying the algebraic machinery to regain the solution previoulsy 
determined by analysis techniques. More than that: The precise form of accessible 
forcing functions is now fully settled, whereas the regularity assumption in ([^ 
was left somewhat vague (a sufficient condition whose necessity was left open). 

3 Two-Point Boundary Problems with One Singularity 

Let us now address the general question of specifying and solving boundary 
problems (as usual: relative to a given fundamental system) that have only one 
singularity. The case of multipliple singularities is left for future investigations. 
Using a scaling transformation (and possibly a reflection), we may thus assume 
the same setting as in Section with the singularity at the origin and the other 
boundary point at 1. 

For the scope of this paper, we shall also restrict ourselves to a certain sub¬ 
class of Stieltjes boundary problems {T,B): First of all, we shall allow only 
local boundary conditions in B. This means multi-point conditions and higher- 
order derivatives (leading to ill-posed boundary problems with distributional 
Green’s functions) are still allowed, but no global parts (integrals); for details 
we refer to [H Def. 1]. The second restriction concerns the differential oper¬ 
ator T G C(x)[9], which we require to be Fuchsian without resonances. The 
latter means the differential equation Tu = 0 is of Fuchsian type (the singu¬ 
larity is regular), and has fundamental solutions x^^(pi{x),..., x^^ with 
each ipi G G‘^{I) having order 0, where Ai,..., A„ are the roots of the indicial 
equation [51 p. 127]. In other words, we do not require logarithms for the solu¬ 
tions. (A sufficient—but not necessary—condition for this is that the roots Ai 
are all distinct and do not differ by integers.) 
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Definition 2. We call {T,B) a boundary problem with mild singularity if T € 
C(a;)[i9] is a nonresonant Fuchsian operator and B a local boundary space. 

For a fixed (T,B), we shall then enlarge the function space T of Section]^ 
by adding ,..., x''" as algebra generators, where each is the fractional 
part of the corresponding indicial root A^. Every element of T is then a sum of 
series J2n>N with /i G {/xi,..., /x„} and N Gh. The integro-differential 
structure on (J^, 9, J) is determined by setting dx^ = XiX^~^, as usual, and by 
using for J the integral as we did also in Section 

As boundary functionals in B, we admit derivatives (0 < ^ < 1, ^ > 0) 
and coefficient functionals {x^^^) {k G Z) whose action is x^ J2n>N “fc- 

For functions / G C^[I) we have of course {x^^^)x^f = /(^^(0)/fc!. Since the 
projectors reg and pp can be expressed in terms of the coefficient functionals, 
we shall henceforth regard the latter just as convenient abbrevations; boundary 
spaces are always written in terms of and {x^) only but of course they can be 
infinite-dimensional. For example, in Sectionj^we had the “regularized boundary 
condition” A(xx) = 0, which is equivalent to Lu = 0 and pp(?x) = 0 and hence 
to {x^) XX = 0 (/c < 0). Its full boundary space is therefore B = [i?, {x^) | /c < 0]. 

The first issue that we must now address is the choice of suitable boundary 
conditions: Unlike in the “smooth case” (without singularities), we may not be 
able to impose n boundary conditions for an rx-th order differential equation. 
The motivating example of Section was chosen to be reasonably similar to 
the smooth case, so the presence of a singularity was only seen in replacing 
the boundary evaluation L by its regularized version A. As explained above, we 
were effectively adding the extra condition pp(xx) = 0 to the standard boundary 
conditions xx(0) = xx(l) = 0. In other cases, this will not do as the following 
simple example shows. 

Example 3. Consider the nonresonant Fuchsian differential equation Tu{x) := 
u” + ^ u' + ^ u = 0. Note that here the indicial equation has the roots Ai = —2 
and A 2 = —1, which differ by the integer 1. Nevertheless, we may take {^, 
as a fundamental system, so T is indeed nonresonant. 

Trying to impose the same (regularized) boundary space B = [pp, L, R] as in 
Example one obtains 

00 00 

nsp = {1 = 0 } 

n——l n——l 

after a short calculation. But this means that the forcing functions / in the 
boundary problem 

u" + ^u' + ^u = f 
pp(xx) = xx(0) = xx(l) = 0 

must satisfy an awkward extra condition (viz. the one on the right-hand side 
of T{B^) above). This is not compensated by the slightly enlarged generality of 
allowing / to have a simple pole at a; = 0. 
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In the present case, we could instead impose initial conditions at 0 so that B = 
[pp, L, LD], In this case one gets T{B^) = C^{I), so there is a unique solution 
for every analytic forcing function. 

For a given nonresonant Fuchsian operator T € C(a:)[9], a better approach 
appears to be the following (we will make this more precise below): 

1. We compute first some boundary functionals /3i,... ,/3„ that ensure a reg¬ 

ular subproblem {T,B) with Bn '■= [Pi,.. ■, fdn]- Adding extra conditions 
(vanishing of all for sufficiently small k) we obtain a boundary 

space B — Bn • such that (T, B) is semi-regular. 

2. If a particular boundary condition is desired, it may be “traded” against one 
of the Pi] if this is not possible, it can be “annexed” to the extra conditions. 
After these amendments, the subproblem {T,Bn) is still regular, and {T,B) 
still semi-regular. 

3. Next we compute the corresponding accessible space T{B^). This space 

might not contain as we saw in Example above when we insisted 

on the conditions u(0) = u(l). 

4. We determine a complement £ of T{B^) as exceptional space in {T,B,£). 

Once these steps are completed, we have a regular generalized boundary 
problem (T, B, £) whose Green’s operator G can be computed much in the same 
way as in Section]^ In detail, we get G = GQ, where G is the Green’s operator 
of the regular subproblem {T,Bn) and Q the projector onto T{B^) along £. As 
we shall see, the operators G and Q can be computed as in the usual setting HU. 
Let us first address Step 1 of the above program. 

Lemma 4. Let T G C(x)[i9] be a nonresonant Fuchsian differential operator of 
order n. Then there exists a fundamental system ui,...,Un € F of T and n 
coefficient functionals Pi := ... ,Pn '■= ordered as ki-\- pi < 

■ ■ • < kn -\- pin so that P{u) G ig a lower unitriangular matrix. 

Proof. We start from an arbitrary fundamental system 

Ui X^ ^ ^ ^l.k^ , ■ • • , ^ ^ ^k'>kn^ 

k'>ki k^kji 

of the Fuchsian operator T, where we take /ii,..., fracational as before and 
we may assume that ai^ki, ■ ■ ■, 0‘n,k„ = 1 so that each fundamental solution Ui 
has order ki. (The order of a series u = x^ 'l2k>N i® defined as the small¬ 
est integer k such that {x^)u 0.) We order the fundamental solutions such 

that ki-\-Pi < ■ ■ ■ < kn-\-pn- We can always achieve strict inequalities as follows. 
If * < n is the first place where ki-\-pi = ki+i -I- pi+i we must also have pi = pi+i 
since 0 < pi,pi+i < 1. Therefore we have ki = ki+i, and we can replace Ui+i 
by Ui+i — Ui and make it monic so as to ensure ki-\- pi < ki+i Pi+i. Repeating 
this process at most n —1 times we obtain ki-\-pi < ■ ■ ■ < kn-\-pn- Choosing now 
the boundary functionals Pi := ..., Pn ■= (a;^"+^'") as in the statement 

of the lemma, we have clearly Pi(ui) = I and Pi(uj) = 0 for j > i as claimed. 
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In particular we see that En := (/3i,... ,/3„)(rti,... ,Un) has unit determi¬ 
nant, so it is regular. Setting Bn := [/3i,...,/?«], we obtain a regular boundary 
problem {T,Bn)- Note that some of the may coincide. For each fj,GM:= 
{/ii,..., /r„} let be the smallest of the ki with fi = fii. We expand the n bound¬ 
ary functionals by suitable curbing constraints to the full boundary space 

B:=Bn + \^^&M,k<k^] (5) 

since the inhomogeneous solutions should be at least as smooth (in the sense of 
pole order) as the homogeneous ones. Note that (T, B) is clearly a semi-regular 
boundary problem. This achieves Step 1 in our program. 

Now for Step 2. Suppose we want to impose a boundary condition /3, assuming 
it is of the type discussed above (composed of derivatives and coefficient 
functionals for fractional parts p, of indicial roots). We must distinguish 

two cases: 

Trading. If the row vector r := /3(ui)i=i^,,._„ S C" is nonzero, we can ex¬ 
press it as a C-linear combination ciri -|- • • • -|- c„r„ of the rows ri,..., 
of En- Let k be the largest index such that 0. Then we may express 
as a C-linear combination of r and the remaining rows Vi {i ^ k), hence 
we may exchange Pk with /3 without destroying the regularity of En = 
(/3l, ■ . ■ , Pn){'U‘i, ■ ■ ■ , Un)’ 

Annexation. Otherwise, we have Ker(T) < /3-*-. Together with Ker(T) + Bn — 
T this implies ([/?] n BnY' = jd-^+Bn = ^ and hence (3 ^ Bn hy the identities 
of [71 App. A]. Furthermore, Lemma 4.14 of jl] yields 

T{Bn + 13^)= T{B^ n/3^) = T{B^) n T(/3^) = T(/?^) 

since we have T{Bn) = E from the regularity of {T,B). But this means that 
adding /3 to B as a new boundary condition necessarily cuts down the space 
of accessible functions unless /3 happens to be in the span of the curbing 
constraints (a:^+^) S B added to Bn in 

In the sense of the above discussion (see Example , the first case signifies a 
“natural” choice of boundary condition while the second case means we insist 
on imposing an extra condition (unless it is a redundant curbing constraint). 
Repeating these steps as the cases may be, we can successively impose any 
(finite) number of given boundary conditions. This completes Step 2. 

For Step 3 we require the computation of the accessible space T{B^). We 
shall now sketch how this can be done algorithmically, starting with a finitary 
description of the admissible space B^ as given in the next proposition. The proof 
is unfortunately somewhat tedious and long-wided but the basic idea is simple 
enough: We substitute a series ansatz into the boundary conditions specified 
in B to determine a number of lowest-order coefficients. The rest is just some 
bureaucracy for making sure that everything works out (most likely this could 
also be done in a more effective way). 
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Proposition 5. Let (T,B) be a semi-regular boundary problem of order n with 
mild singularity such that (T, Bn) is a regular subproblem. Let M = {/ri,..., /r„} 
be the fractional parts of the indicial roots for T. Then we have a direct decom¬ 
position of the admissible space B^ = with components 

A^, = [p^,i{x),...,p^iJ^x)]+P^{C^{I)^) ( 6 ) 

Here p^i(x),... ,p^i^{x) £ C[a:, are linearly independent Laurent polynomials 
and the linear operators —>■ C‘^{L) are defined by 

Pf^{b{x)) = x^i^b^ix) -f 


with £ Z and Laurent polynomials q^u^j{x) £ C[a;, almost all of which 
vanish over the summation range 0 < ^ < 1 and j > 0. 

Proof. Any element u € P may be written in the form 

OO 

fiGM k—k(fj,) 

where the k{fi) £ Z are a priori arbitrary. However, because of the curbing 
constraints of P in ([^ we may assume fc(/i) = for u £ B^. Let us write B+ 
for the extra conditions that were “annexed” at Step 2 (if none were added 
then clearly B+ = O); thus P is a direct sum of Bn and B+ and the curbing 
constraints. If denotes the largest integer such that occurs in any 

condition of Bn + B+ and if we set := -|- 1 -f /i, we can also write u £ B-^ as 

Sfj . 

u = Y afj,kx’^+^-\- Y (7) 

fiGM k—k/j, /tGM 

where the 6^ (a:) £ C‘^{I) are convergent power series. Next we impose the con¬ 
ditions f3 G Bn -\- B+ on u. But each /3 is a C-linear combintions of coefficient 
functionals with u >->■ and derivatives (0 < ^ < 1,^ > 0) whose 

action on u yields 


E E +E E H) ( 0 , 

fi k (I j=0 ^ ^ 

Hence j3{u) yields a C-linear combination of the a^k and the b^J\^). We compile 
the coefficients £ C into a column vector a £ C'^ of size S := + 

consisting of |M| contiguous blocks of varying size s^ —A:^-|-l. Putting the (finite) 
set S of evaluation points occurring in Bn +P+ into ascending order, we compile 
also the derivatives 6^(^), 6(j(^),... into a column vector b £ of size T := 
consisting of \M\ contiguous blocks, each holding derivatives b^J'^ (f) for 
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a fixed ^ G M and certain ^ G ^ and j > 0. Assuming Bn+B+ has dimension R, 
the boundary conditions /3(u) =0 for /3 G can be written as 

Aa = Bb (8) 

for suitable matrices A G and B G that can be computed from the 

boundary functionals f3 G Bn + B+. 

Regarding the right-hand side as given, let us now put Q into row echelon 
form (retaining the same letters for simplicity). If the resulting system contains 
any rows that are zero on the left but nonzero on the right, this signals constraints 
amongst the b^{x) whose treatment is postponed until later. For the moment we 
discard such rows as well rows that are zero on both sides. Let U < R he the 
number of the remaining rows andpi,... ,pu the pivot positions and V := S — U 
the number of free parameters. Then we can solve ([^ in the form a = a + A-C^. 
Here a G C'® is the vector with entry {Bb)j in tow pj and zero otherwise, and A G 
(j^Sxv consists of the non-pivot columns of the corresponding padded matrix 
(its j-th row is the i-th row of A for pivot indices j = Pi, and —ej for non¬ 
pivot indices j). Writing the free parameters as D = (ui,..., uy) G C^, we may 
substitute this solution a into the ansatz Q to obtain 

■u = ^ ^ -k ^ 6^(x), (9) 

fieM k=kij, fieM 

where the a^kiv,b) are C-linear combinations of the free parameters v and the 
derivatives b^J\^) comprising b. Note that we may regard the a^fc(...) as row 
vectors in whose entries are computed from A and B. 

We turn now to the issue of constraints amongst the b^{x), embodied in those 
rows of the reduced row echelon form of Q that are zero on the left-hand side 
but nonzero on the right-hand side. We put the corresponding block of B into 
reduced row echelon form (again retaining the same letters for simplicity). Let p,' 
the smallest u occurring in any such row. Then each of the rows containing a 
derivative b^^} (x) provides a constraint of the form 

H H ’ (0 = 0, (10) 

ies j {gh j 

where the cj^-, G C are determined by B, and with finite sums over j. Let the 
number of such constraints be X. Collecting the into a vector b^i G 

of size Y ■.= we can write the constraints ( |10[ ) as matrix equation Cb^i = 
d where the coefficient matrix C G is determined by the c^j of each 

constraint ( |l0| ), and the right-hand side d G by the corresponding and 
the ^ (^) for p > /i', which for the moment we regard as known. Note that X < 
Y since C is in row echelon form. Let its pivots be in the positions qi,... ,qx 
and set Z := Y — X. Then we can write the solution as b^i = 6^/ -k (7 • . 

Here 6^/ G is the vector with entry dj in row qj and zero otherwise, while C G 
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consists of the non-pivot columns of the padded matrix (its j-th row is 
the i-th row of C for pivot indices j = qi, and —Cj for non-pivot indices j). 
Writing w = (tci ,... ,wz) & for the corresponding free parameters, we may 
thus view (10) as providing Y constraints 




( 11 ) 


where the b^>^j(w,b^) are C-linear combinations of the free parameters w and 
certain derivatives b^ii\^) that we have collected into a vector 6_|_ S . 

Again we may view b^i^j [...) as a row vector in whose entries can 

ultimately be computed from A and B. 


We regard now (111 as determining equations for fixing Y of the coefficients 
of bfj,i {x) = b^'kx'^- Indeed, if j is the highest derivative order occurring in (111 

we may split according to b^i{x) = J2k<j b^/kX^ + x^'^^b^i{x), with an arbitrary 
power series b^i(x), and then substitute this into 0 to obtain 


Vj = - 0 


i=0 


3 + 1 
i + 1 


i+i rb 




( 12 ) 


for fixing one coefficient of 6^/(a;). Now we substitute this bj back into the above 
ansatz 6^/ (a;) = J2k<i b^ikX^ Yx^^^b^i (a;), and we repeat the whole process for all 
other constraints ( |11[ ), each time determining the lowest unknown coefficient bk 
of b^i{x). Of course, it may be necessary to expand the splitting to extract a 
larger polynomial part (if all the coefficients in the current polynomial part are 
determined). Eventually, we end up with 


bii'{x)= ^ b^'k{w,h+)x'" + x'^t^'^'^bf,,{x), (13) 

k<m^r 


with some break-off index m^i that is specific to 5^/(x), and with 6+ enlarged 
to comprise also the derivatives b)'^) {C) that were needed in determining the 
coefficients (12|. Substituting (13) into ([^, we see that the b^ik(w,b+) may 


be combined with the a^ik{v,b) if we adjoin the parameters w to the parame¬ 
ters D; then we can also rename the series bfi'{x) back to bf^i{x). Of course new 


terms a^ik{v, b) x^'^^ may be created in the ansatz ([^, and its polynomial part 
may be expanded—but its overall form is not altered. 

We have now eliminated those constraints amongst the {x) occurring in ([^ 
that involve b^i, where fj,' was chosen minimal. Hence we are only left with 
constraints amongst the b^{x) with ^ . Repeating the elimination process 

a finite number of times (at most \M\ eliminations are necessary), we obtain 
the generic form of u G as given in (§, where the b^(x) are now arbitrary 
(convergent) power series. Since the terms with distinct factors x^ are clearly in 
distinct direct sum components x^A^^ the latter may now be described by 


Afi 


^ a^fc(h, b) x’^ + x^i^ b^{x) 

k=kfj. 


VGC^, \ 

b{x) G C'‘^(/)^ j’ 
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where we have set := +1. Splitting the coefficients h) into a C-linear 

combination of the free parameters (ui, ... ,vv) and a C-linear combination of 
the derivatives comprising h, we collect terms whose coefficient is a specific 

parameter Vi {i = 1,... ,V) or a specific derivative hv\^) 0). 

This leads to 

( y 

^ i—1 h>^M 

and hence to by extracting a C-basis for each [p^i(a;),... ,Pf_tviA]- 

We have now established Step 3 of our program since the accessible space T{B^) 
can be specified by applying T G C[x,^ to the generic functions ^ of the ad¬ 
missible space B-^. Our next goal is to find a projector Q onto T{B^), which 
then gives the exceptional space £ := Ker((5) = Im(l — Q) required for Step 4. 
This will be easy once we have a corresponding projector onto In fact, the 
operator in Proposition is not quite a projector (for one thing, it is not 
even an endomorphism), but in a sense it is not far away from being one. For 
seeing this, note that we have 

J-= 0 cr^C((cr)), (14) 

mgm 



where C((a:)) denotes the field of Laurent series (converging in the punctured 
unit disk). The direct decomposition (14) just reflects the fact that the for 
distinct p G M are linearly independent. Let us write {p) : P ^ F ioi the indicial 
projector onto the component x^C((x)) of (14). In other words, (p) extracts all 
terms of the form (k G Z) from a series in P. Note that combinations 
like P = e^D^{p) provide linear functionals P G P* for extracting derivatives of 
the /r-component of a given series in P. 

For writing the projector corresponding to Proposition!^ let us also introduce 
the auxiliary operator x^: C((x)) —)■ x'^C((x)) < P and its inverse x~^. Then 
we shall see that the required projector is essentially a “twisted” version of two 
kinds of projector: one for splitting off the polynomial part of the occurring 
series, and one for imposing the derivative terms. For convenience, we shall use 
orthogonal projectors in C((x)), where the underlying inner product is defined 
by (x^|x^) = 5ki for all k,l G Z. Such projectors are always straightforward to 
compute (using linear algebra on complex matrices). 


Corollary 1. Using the same notation as in Proposition^ let i?^, : C(x)) —)■ 

C(x) be the orthogonal projectors onto [p^i(x),... ^p^ijpxp and onto xP C((x)), 
respectively. Writing R'^ := x^R^x~^{p) and S'^ := x^Sfj_x~^(p) for their twisted 
analogs, we define the linear operator P: P ^ P by 


^ - X {k + X X 

AtGM 


is a projector onto . 


(15) 
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Proof. The action of P on a general series (which we may split at = J/x ~ 1 
in its /r-components) 


u{x) = f^{x) =Y^ + Y bf,{x) e P 

fL£M k—k^ 

is 

Pu{x) = Y Y ft^kx'^) 

fiGM ^ k—kfj^ ^ 

where we have set bi_i{x) := x~P f^{x) € C‘^{I). Extracting the x^ component 
of Pu yields 

Rri^Y ft^kX^) + 2^^"“ H XI e Af, (16) 

k=kf_, ueM ^,j 


as one sees by comparing with the last displayed equation in the proof of Propo¬ 
sition]^ Hence we may conclude that Im(P) < B^. 

It remains to prove that Pu = u for u G since this implies Im(P) > B^ 
and P^ = P so that P is indeed a projector onto B^ as claimed in the corollary. 
So assume u{x) G B^ is arbitrary, and split it as above. The orders of the se¬ 
ries ffi are given by the curbing constraints. Since u{x) already satisfies all bound¬ 
ary conditions in -|- B+ < B, the reduced row echelon form of Q will contain 
only zero rows. The original ansatz (J^ is then left intact; no expansion of the 
polynomial part is necessary and no b are involved in its coefficients. Therefore 
the original coefficients a^k in 0 coincide with the coeffients a^k{v,b) = a^k{v) 
in Q, which constitute the polynomials of := [p^i{x),... ,p^i^{x)]. Now let 
us consider (16). Since the series ffj,kX^ is thus already in = Im(i?^), we 
may omit the action of and since a^k{v, b) = o^fe('()) the triple sum in ( [l^ is 
zero. But then Pu{x) becomes identical with u{x) as was claimed. 


For accomplishing Step 4 of our program, it only remains to determine a 
projector Q onto the accessible space T{B-^) from the projector P onto the 
admissible space B^ provided in Corollory]^ This can be done easily since Q is 
essentially a conjugate of P except that we use a fundamental right inverse T^, 
for want of a proper inverse. (The formula for in [TTl Prop. 23] and [TU 
Thm. 20] may be used but recall that in our case J = so that e = Ei.) 

Proposition 1. Using the same notation as in Proposition^ the operator 


Q := TPT^ : P^P 


is a projector onto the accessible space T{B^). 
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Proof. Let us first check that Q is a projector. Writing U := 1 — T^T for the 
projector onto Ker T along [e, eU, ..., eZ 3"“^] and using = P, we will indeed 
get 

q 2 ^ yp2y0 _ XPUPT^ = TPT^ = Q 

provided we can ascertain that KerT < KerP =: C since then PU = 0. We 
know that + C = P since P is a projector onto along C. On the other 
hand, we have also Ker(T) + = P since (T,Bn) is regular. Intersecting C in 

the former decomposition with the latter yields 

B^ + {CnB:^) + {CnKeTT)=P. (17) 


But Bn < B implies B^ < ^n’’ intersecting the decomposition B^ C = P 
with Bn leads to B-^ + (C n B^) = Bn - Using the other decompostion Ker(T) + 
Bn = P one more time we obtain 

P^ + (CnP^) + KerT = P. (18) 


Comparing 0 and ( |I^ , we can apply the well-known rule [H (2.6)] to obtain 
the identity C C KerP = KerP and hence the required inclusion KerP < C. 

It remains to prove that Im((5) = T{B^). The inclusion from left to right is 
obivous, so assume f = Tu with u € B^. Then f = u—Uu and hence PT^ f = 
Pu — PUu = u because P projects onto B^ and PU = 0 from the above. But 
then we have also Qf = TPT^ f = Tu = f and in particular / G Im((5). 


We have now sketched how to carry out the four main steps in our program 
aimed at the algorithmic treatment of finding/imposing “good” boundary con¬ 
ditions on a Fuchsian differential equation with one (mild) singularity. At the 
moment we do not have a full implementation of the underlying algorithms in 
TH30REMV (or any other system). However, we have implemented a prototype 
version of some portion of this theory. We shall demonstrate some of its features 
by with example from engineering mechanics. 


4 Application to Functionally Graded Kirchhoff Plates 

Circular plates play an important role for many application areas in engineering 
mechanics and mathematical physics. If the plates are thin (the ratio of thickness 
to diameter is small enough), one may employ the well-known Kirchhoff-Love 
plate theory [^, whose mathematical description is essentially two-dimensional 
(via a linear second-order partial differential equation in two independent vari¬ 
ables). We will furthermore restrict ourselves to circular Kirchhoff plates so as 
to have a one-dimensional mathematical model, via a linear ordinary differential 
equation of second order. 

However, we shall not assume homogeneous plates. Indeed, the precise man¬ 
ufacture of functionally graded materials is an important branch in engineering 
mechanics. In the case of Kirchhoff plates, the functional grading is essentially 
the variable thickness t = t{r) or variable bending rigidity D = D{r) of the plate 
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along its radial profile. (We write r for the radius variable ranging between zero 
at the center and the outer radius r = b.) 

Let w = w(r) be the displacement of the plate as a function of its radius. 
This is the quantity that we try to determine. It induces the radial and tangential 
moments given, respectively by 

Mrir) = -D(r) {w"{r) + ^w'{r)), 

Mg{r) = -D{r) w'{r) + vw''(r)), 


where v is the Poisson’s ratio of the plate (which we assume constant). For 
typical materials, v may be taken as g, jj | or even 0. A reasonable constitutive 
law for the bending rigidity is 


D{r) 


E(r)t{rf 


(19) 


where E = E{r) is the variable Young’s modulus of the plate. 
The equilibrium equation can then be written as 


dMr Mr - Me 
dr r 


Qrt 


where Qr = Qr('c) is the cumulative load 


( 20 ) 


1 r 1 r 

Qr{r) = -/ q{r)2Trrdr = -/ q{r)rdr 

27rr r Jq 


induced by a certain loading q = q{r) that may be thought to describe the weight 
(or other forces) acting in each ring [r, r -|- dr]. 

For the calculational treatment of (20) it it useful to introduce the func¬ 
tion (fi := —w'(r), which represents the (negative) slope of the plate profile. In 
terms of q), the equilibrium equation is given by 


A typical example of a useful thickness grading is the linear ansatz t = tQ{l — ^), 
cut off beyond some a < b close to b; this describes a radially symmetric pointed 
plate with straight edges (more or less a very flat cone). Suppressing the cut-off 
for the moment and changing the independent variable r to p = r/6, we have thus 
thickness t{p) = to ■ {1 — p) and from (19) bending rigidity D{p) = Dq ■ {1 — pY 
with Dq := Eto/12{l — v^). The equilibrium equation becomes now 



1 \ y’(p) 

I- p) p 


Qrip) b"^ 
Do{1-pY^ 


where we have set v = ^ for simplicity. Note that the right-hand side of this equa¬ 
tion, which we shall designate by f{p), is not a fixed function of p but depends 
on our choice of the loading. Hence we consider f{x) as a forcing function. 
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For the boundary conditions (for once we use this word in its original literal 
sense!) we shall use u;'(0) = w'{a) = 0, which translates to </9'(0) = = 0 

in the ip = ip(p) formulation, with the abbreviation /3 := a/b < 1. Physically 
speaking, this corresponds to a plate that is clamped in the center and left free 
at its periphery (this comes from translating suitable boundary conditions for 
the displacement w = w{r) and taking the appropriate limits). In summary, we 
have the boundary problem 


t="(rt + (G A) 

(p(0) = (p(/3) = 0, 

P'(p) - 


^=f(p), 


( 21 ) 


which is indeed of the type discussed in Section by a simple scaling from / = 
[0,1] to [0,/?]. Its treatment in the GreenGroebner package proceeds as follows 
[we apologize for the lousy graphics rendering—we plan to fix this as soon as 
possible]: 


soli: BPSolvej 

♦ " [p] . (2 - -t-U' (Pi - ) ♦(pl = £[pl 

\p l-pl \p’ (I-P)P/ 

♦ [L] =*(R] irO 


fi, Paraas -• (B}< PoetProc -» PullSiinplify| 

1 1 

6 (-11 

3)^ P 





(pb-. 

t2pl [ 

-t(3.2plp 

-t3(. 

-2plp' 

t3p'(-3t2p) {fIdJt - 


J 

{-ltd' 

2" (-ltd' 


• i-ltd' iM-itd' 

3 


-{blflPi t3(3-2p)p' 

" -til 


;'fidii{ tp'i-3t2pi ['-{'([{M? - 



i-itd' 

• (-ltd' 

(-ltd' 

15 

r- 

, fP 1 

-9 -: 

, fP 1 

.2 -; 


1 11 
t(3-2plp' -f'fldtif (-3-3-It 



■'« (-ltd' 


k I-ltd' p-’ -3t2F 

(3 

2PIP^ 

' -{bifidf (Itll 

--|t(3-2plp' 1 

—- 

bldiff 1-15-It 



• (-ltd' 

02-3 + 2^ Jo 

(-ltd' 

.3 + 2^ 

(3 

2 pip’ 

' -(-2- 

—l + {3-2p)p^ [^- 

— dfifipf 13-)i 



" (-ltd' p 

-3 + 2p Jo {. 

td' 

^2 -3 + 2P 


The output is the Green’s operator G of (21) applied to a generic forcing 


function /(p), giving the solution ip(p) = Gf (p) as an integral 


1 


9{p,0fi0dC 


in terms of the Green’s function 5 (p, ^); the latter can also be retrieved explicitly 
if this is desired (note that the expression is clipped on the right-hand side so 
the two case conditions ^ < p and p < ^ labelling the two lines are not visible): 


, ♦, p, 0, s, Panis ^ (f}, PoatProc ^ FullSinplify, 


Output -t GreensFunction 


pf: BPSolve1 

F'[pl*(---U'[pl-(^ t—)^(pl---f[pl 

\p l-p/ \p‘ |l-p)p/ 


♦ [L! = *[I!1 = 0 
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Let us now choose a constant loading q{p) = qq. Taking the cut-off into 
account, this leads to the forcing function f{p) = In this case, we can 

compute the solution ip{p) of (21) explicitly. For definiteness, let us choose a 
cut-off at /? = 0.9. In that case, we obtain 


ip(p) = (^- 2790p3 -f 1944 log(9/p - 9) + 2000 log(l - p) 

- 2000p3 log(10 - lOp) -f 4671p2 - 2916p2 log(9/p - 9) - 3000p^ log(l - p) 
-f 3000p^ log(10 - lOp) - 1944p -f 972 log(l - p)) j (2916(p - l)^p), 

and the corresponding displacement w{p) = —Jg (p(p) dp is given by 

w{p) = ( - 972(p - 1) Li2 ((1 - p)-i) -f 972(p - 1) Li2(l - p) - 2790p2 
-f 1944p2 log(9) -f 3944p2 log(l - p) - 2000p2 log(-10(p - 1)) 

- 1944p2 log(p) -f 2853p -f 500plog^(10) -f 14plog^(l - p) 

- 500plog^(-10(p- 1)) - 141og2(l - p) -f 5001og^(-10(p - 1)) 

- 972plog(9) - 1500plog(100) -f 486plog(81) log(l - p) - 7825plog(l - p) 
-I- 4000plog(-10(p - 1)) -f 972plog(l - p) log(p) -f 972plog(p) 

- 1944 log (9/p — 9) — 4861og(81) log(l — p) + 6825 log(l — p) 

- 30001og(-10(p - 1)) - 9721og(l - p) log(p) - 19441og(p) - 5001og^(10) 
-f 19441og(9) -f 15001og(100)) j 2916(1 - p) 


We have displayed the graphs of these solutions (p(p) and rc(p) below. 
9(P) W(P) 




Of course one could use different functional gradings and/or loading func¬ 
tions, and the integrals would not always come out in closed form. In this case one 
can resort to numerical integration (which is also supported by Mathematica). 
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